function [tOut uOut] = hydro_solve(params, uIn)

  time_steps = params.t_steps;


  global state_t_prev;
  global state_n;
  global prev_t;
  global fps; 
  state_t_prev = 0;
  state_n = 0;
  prev_t = 0;
  fps = 0;
 
  global xt
  xt = [];
  clear xt;
  global xt;
  xt = xtimer();

  % ode45 options
  
  if params.InitialStep ~= []
    opts = odeset('RelTol', params.RelTol, 'AbsTol', params.AbsTol, 'OutputFcn', @(t,u,flag) hydro_plot(t, u, flag, params), 'InitialStep', params.InitialStep);
  else
    opts = odeset('RelTol', params.RelTol, 'AbsTol', params.AbsTol, 'OutputFcn', @(t,u,flag) hydro_plot(t, u, flag, params));
  end

  %  2 -- 2nd order hydrodynamics, see mmatica file
  % -2 -- 2nd order hydro with terms O(ux,uy)^2 set to zero before inversion. the nr code is not as optimized as rel 2nd order and might actually run slower
  %  0 -- 0th order hydrodynamics
  switch params.hydroOrder
    case 2
      eom = @hydro_eom2;
    case -2
      eom = @hydro_eom2nr;
    case 0
      eom = @hydro_eom0;
  end

  hydro_begin=tic;
  % run the ode solver
  %[tOut, uOut] = ode45(@(t,u) eom(t,u,params), time_steps, uIn, opts);

  if length(time_steps) == 2
    sol = ode45(@(t,u) eom(t,u,params), time_steps, uIn, opts);
    tOut = sol
    uOut = 0;
  else
    [tOut uOut] = ode45(@(t,u) eom(t,u,params), time_steps, uIn, opts);
  end

  disp(sprintf('computed hydro simulation in %fsec', toc(hydro_begin)));

  xt.stats();
  return


  % old code, before changing to ode sol

  % % ****************************
  % % compute 1st time derivatives
  % % ****************************

  % % first derivatives computed using the hydro eom code itself. 
  % % pity ode45 doesn't return that by default.
  % % cost to compute the 1st derivatives at pre-set time points is insignificant in comparison with the original solving of the ode

  % tic;
  % Nt = params.Nx*params.Ny;
  % udtOut   = zeros(params.Tsteps, Nt.*5);

  % % first derivative from eom
  % for it=1:params.Tsteps
  %   t = tOut(it);
  %   u = uOut(it, :);
  %   
  %   udtOut(it,:) = eom(t, u, params);
  % end

  % % ****************************
  % %  get rid of Pi and B - no longer needed
  % % ****************************

  % uOut  (:,3*Nt+1 : 5*Nt) = [];
  % udtOut(:,3*Nt+1 : 5*Nt) = [];


  % % ****************************
  % % compute 2nd time derivatives
  % % ****************************

  % % second derivative computed using 1st derivate finite-difference scheme

  % derorder = params.timeOrder;
  % derker   = d1ker(params.timeOrder);
  % udtdtOut = conv2(udtOut, derker, 'same')./params.Dt;

  % % delete derivative padding
  % deli = [1:(derorder/2) (length(tOut)-derorder/2+1):length(tOut)];
  % tOut    (deli, :) = [];
  % uOut    (deli, :) = [];
  % udtOut  (deli, :) = [];
  % udtdtOut(deli, :) = [];

  % elapsed = toc;
  % disp(sprintf('computed higher hydro time derivatives in %fsec', elapsed));

  % % adjust begin/end times
  % params.t_beg = params.t_beg + derorder./2.*params.Dt;
  % params.t_end = params.t_end - derorder./2.*params.Dt;
  % params.Tsteps = params.Tsteps - derorder;


  % % ****************************
  % % compute ut and derivatives
  % % ****************************

  % ux = uOut(:, 1*Nt+1 : 2*Nt);
  % uy = uOut(:, 2*Nt+1 : 3*Nt);

  % uxdt = udtOut(:, 1*Nt+1 : 2*Nt);
  % uydt = udtOut(:, 2*Nt+1 : 3*Nt);

  % uxdtdt = udtdtOut(:, 1*Nt+1 : 2*Nt);
  % uydtdt = udtdtOut(:, 2*Nt+1 : 3*Nt);

  % ut = sqrt(1+ux.^2+uy.^2);
  % utdt = (uxdt.*ux + uydt.*uy)./ut;
  % utdtdt = (1+ux.^2+uy.^2).^(-3/2).*(uxdt.^2.*(1+uy.^2)+(-2).*ux.*uxdt.*uy.*uydt+(1+ux.^2).*uydt.^2+(1+ux.^2+uy.^2).*(ux.*uxdtdt+uy.*uydtdt));

  % uOut    (:, 3*Nt+1 : 4*Nt) = ut;
  % udtOut  (:, 3*Nt+1 : 4*Nt) = utdt;
  % udtdtOut(:, 3*Nt+1 : 4*Nt) = utdtdt;



end


function ker = d1ker(order)
  switch (order)
    case 2
      ker = [-1/2; 0; 1/2];
    case 4
      ker = [1/12; -2/3; 0; 2/3; -1/12];
    case 6
      ker = [-1/60; 3/20; -3/4; 0; 3/4; -3/20; 1/60];
    case 8
      ker = [1/280; -4/105; 1/5; -4/5; 0; 4/5; -1/5; 4/105; -1/280];
  end
  ker = flipud(ker);
end


function ker = d2ker(order)
  switch (order)
    case 2
      ker = [1; -2; 1];
    case 4
      ker = [-1/12; 4/3; -5/2; 4/3; -1/12];
    case 6
      ker = [1/90; -3/20; 3/2; -49/18; 3/2; -3/20; 1/90];
    case 8
      ker = [-1/560; 8/315; -1/5; 8/5; -205/72; 8/5; -1/5; 8/315; -1/560];
  end
end


